/*=============================================================================
*
*   Copyright (C) 2020 All rights reserved.
*
*   Filename:		model.hpp
*
*   Author: Wang Zhecheng - wangzhecheng@yeah.net
*
*   Date: 2020-02-16 22:55
*
*   Last Editors: Wang Zhecheng - wangzhecheng@yeah.net
*
*   Last modified:	2020-03-30 21:12
*
*   Description:
*
=============================================================================*/
#ifndef IMODEL_H
#define IMODEL_H

#include "../arithmetic/registrator.hpp"
#include "../geometry/mesh/build_mesh.hpp"
#include "../geometry/registrator.hpp"
#include "../io/data_parser.hpp"
#include "../material/registrator.hpp"
#include "../util/observer_interface.hpp"
#include "../xfem/build_crack.hpp"
#include "../xfem/propagate.hpp"
#include "boundary.hpp"
#include "dof_handler.hpp"
#include "solve/solve_time_step.hpp"
#include "time_phase.hpp"
#include "time_step.hpp"

namespace geoxfem {
namespace model {

template <typename method_policy,    // fem or sth else
          typename couple_policy,    // solid mechanism or Biot of sth else
          typename dimension_type,   // plane stress, plane strain or 3-D
          typename equation_policy,  // static or dynamic
          typename solve_policy,     // implicit or explicit
          typename load_policy,      // incremental or absolute load policy
          typename degree>           // Shape function degree (Linear)
struct Model {
  using MethodPolicy   = method_policy;
  using CouplePolicy   = couple_policy;
  using DimensionType  = dimension_type;
  using EquationPolicy = equation_policy;
  using SolvePolicy    = solve_policy;
  using LoadPolicy     = load_policy;
  using Degree         = degree;
  using SELF           = Model<MethodPolicy,
                     CouplePolicy,
                     DimensionType,
                     EquationPolicy,
                     SolvePolicy,
                     LoadPolicy,
                     degree>;

  using ptree        = boost::property_tree::ptree;
  using Reader       = io::DataParser<io::json>;
  using point_type   = typename traits::point_type<DimensionType>::type;
  using Shape        = geometry::ShapeInterface<DimensionType, degree>;
  using DOFHandler   = DofHandler;
  using Material     = material::MaterialInterface<DimensionType>;
  using Mesh         = geometry::Mesh<DimensionType, degree>;
  using DirichletBC  = Boundary<dirichlet>;
  using NeumannBC    = Boundary<neumann>;
  using TimeFunc     = arithmetic::Time_Function_Interface;
  using TStep        = TimeStep<SolvePolicy>;
  using Observer     = util::ObserverInterface<SELF>;
  using ObserverList = std::list<Observer*>;

  using PointList = std::unordered_map<size_t, point_type>;
  using ShapeSpList =
    std::unordered_set<std::unique_ptr<Shape>, shapeHasher, shapeCmp>;
  using MaterialList = std::unordered_map<size_t, std::unique_ptr<Material>>;
  using Crack        = xfem::Crack<DimensionType, degree>;
  using CrackList    = std::vector<Crack>;

  using Vertex  = geometry::Vertex<DimensionType, degree>;
  using Segment = geometry::Segment<DimensionType, degree>;
  using Polygon = geometry::PolygonInterface<DimensionType, degree>;

  static constexpr size_t dimension = traits::dimension<DimensionType>::value;

  Model()
      : m_time{0},
        m_iter{0},
        m_list_point{},
        m_list_shapes{4},
        m_list_material{},
        m_list_crack{},
        m_dirichlet{},
        m_neumann{},
        m_mesh{},
        m_dof_handler{},
        m_list_phase{},
        m_max_pid{0},
        m_reader{} {}
  ~Model() = default;
  /**
   * @brief read
   *
   * @param filename
   */
  void init(Reader& reader) { /*{{{*/
    m_reader = &reader;
  } /*}}}*/
  /**
   * @brief build this Problem
   */
  void buildProblem() { /*{{{*/
    GEOXFEM_LOG_SEV(info) << "Building problem...";
    /* 0. register everything */
    register_type();
    /* 1. set up constances */
    setupParams();
    /* 2. build material list */
    buildMaterial();
    /* 3. build Phases */
    buildTimeSequence();
    /* 4. create mesh object and populate mesh by xml resource data */
    buildMesh();
    /* 5. create cracks */
    buildInitCracks();
    /* 6. output/plot initial mesh */
    notify();
  } /*}}}*/
  /**
   * @brief solve this Problem
   */
  void solveProblem();
  /**
   * @brief register_type
   */
  static void register_type() { /*{{{*/
    util::registrator<Shape>::apply();
    util::registrator<Material>::apply();
    util::registrator<TimeFunc>::apply();
  } /*}}}*/
  /*direct member access{{{*/

  /**
   * @brief get data reader
   *
   * @return
   */
  Reader* getReader() {
    return m_reader;
  }
  /**
   * @brief getMesh
   *
   * @return
   */
  Mesh* getMesh() {
    return &m_mesh;
  }
  /**
   * @brief getMesh
   *
   * @return
   */
  const Mesh* getMesh() const {
    return &m_mesh;
  }
  /**
   * @brief getDofHandler
   *
   * @return
   */
  DOFHandler* getDofHandler() {
    return &m_dof_handler;
  }
  /**
   * @brief getDirichletBoundary
   *
   * @return
   */
  DirichletBC* getDirichletBoundary() {
    return &m_dirichlet;
  }
  /**
   * @brief getNeumannBoundary
   *
   * @return
   */
  NeumannBC* getNeumannBoundary() {
    return &m_neumann;
  }
  /**
   * @brief get pointer to #n material
   *
   * @param n
   * @return
   */
  const Material* getMaterial(const size_t& n) const {
    return m_list_material.at(n).get();
  }
  /**
   * @brief getTime
   *
   * @return
   */
  _real getTime() const {
    return m_time;
  }
  size_t getIter() const {
    return m_iter;
  }
  void setIter(const size_t& value) {
    m_iter = value;
  }
  /*}}}*/
  /*point manuvar{{{*/

  /**
   * @brief append
   *
   * @param point
   *
   * @return
   */
  point_type append(point_type& point) {
    int id = point.getIndex();
    if(id == -1) {
      id = m_max_pid;
      m_max_pid++;
    }
    else if(id >= m_max_pid)
      m_max_pid = id + 1;

    point.setIndex(id);
    auto [it, success] = m_list_point.emplace(id, point.deepClone());
    return it->second.shallowClone();
  }
  /**
   * @brief append a smart pointer of Point
   *
   * @param spPoint
   *
   * @return
   */
  point_type append(point_type&& point) {
    int id = point.getIndex();
    if(id == -1) {
      id = m_max_pid;
      m_max_pid++;
    }
    else if(id >= m_max_pid)
      m_max_pid = id + 1;

    point.setIndex(id);
    auto [it, success] = m_list_point.emplace(id, std::move(point.deepClone()));
    return it->second.shallowClone();
  }
  /**
   * @brief getPoint
   *
   * @param n
   *
   * @return
   */
  point_type getPoint(const size_t& n) {
    return m_list_point.at(n).shallowClone();
  }
  /*}}}*/
  /*shape manuvar{{{*/

  /**
   * @brief append a smart pointer of shape
   *
   * @param spShape
   *
   * @return
   */
  Shape* append(std::unique_ptr<Shape>& spShape, Mesh* mesh = nullptr) {
    Shape* shp      = spShape.get();
    size_t topo_dim = shp->getTopologyDimension();
    int    id       = shp->getIndex();
    id              = (id == -1) ? m_list_shapes.at(topo_dim).size() : id;
    shp->setIndex(id);
    std::vector<Shape*> vlist;
    auto [it, success] = m_list_shapes.at(topo_dim).emplace(std::move(spShape));
    if(success) {
      /* 1. append facets of shape */
      auto facets = shp->createFacet();
      for(auto& facet : facets) {
        Shape* v = append(facet, mesh);
        vlist.emplace_back(v);
      }
      /* 2. link two vertices of segment */
      if(topo_dim == 1) {
        vlist.front()->link(*vlist.back());
        vlist.back()->link(*vlist.front());
      }
      /* 3. build shape connectivity and boundary status */
      buildConnectivity(*shp);
      /* 4. build standard gauss points */
      model::buildGaussPoints(*shp);
      /* 5. add shape to given mesh */
      if(mesh)
        mesh->append(shp);
    }
    return it->get();
  }
  /**
   * @brief getNumberOfShapes
   *
   * @param dim
   *
   * @return
   */
  size_t getNumberOfShapes(const size_t& dim) const {
    return m_list_shapes.at(dim).size();
  }
  /**
   * @brief search shape by dimension and id. terrible way to access shape. Only
   * used in building node-element connectivity.
   *
   * @param dim
   * @param id
   *
   * @return
   */
  Shape* getShape(const size_t& dim, const size_t& id) const {
    auto& list    = m_list_shapes.at(dim);
    auto  find_id = [id](const std::unique_ptr<Shape>& shp) -> bool {
      return shp->getIndex() == id;
    };
    auto it = std::find_if(list.begin(), list.end(), find_id);
    return it->get();
  }
  /**
   * @brief buildConnectivity
   *
   * @param shp
   */
  void buildConnectivity(Shape& shp) const {
    size_t topo_dim = shp.getTopologyDimension();
    /* 1. build connectivity to vertices; */
    if(!shp.ifVertex())
      for(auto& p : pointsIn(shp)) {
        size_t vid = p.getIndex();
        Shape* v   = getShape(0, vid);
        shp.link(v);
        v->link(shp);
      }
    /* 2. build connectivity to rest shapes */
    for(size_t i = 1; i != 4; ++i)
      for(auto& it : shapesIn(*this, i)) {
        auto&              shp2     = *it;
        geometry::topology relation = shp2.getTopologyWith(shp);
        if(relation == geometry::topology::child
           || relation == geometry::topology::parent
           || relation == geometry::topology::adjacent) {
          shp.link(shp2);
          shp2.link(shp);
        }
      }
    /* 3. update boundary status */
    if(topo_dim > 0)
      for(auto& facet : neighborsIn(shp, topo_dim - 1))
        facet->updateBoundaryStatus();
    return;
  }
  /*}}}*/
  /*crack manuvar{{{*/
  Crack* append(Crack& crack) {
    m_list_crack.emplace_back(crack);
    Crack* p_crack = &m_list_crack.back();
    p_crack->init();
    return p_crack;
  }
  Crack* append(Crack&& crack) {
    m_list_crack.emplace_back(std::move(crack));
    Crack* p_crack = &m_list_crack.back();
    p_crack->init();
    return p_crack;
  }
  Crack& getCrack(const size_t& n) {
    return m_list_crack.at(n);
  }
  const Crack& getCrack(const size_t& n) const {
    return m_list_crack.at(n);
  }
  CrackList& getCrack() {
    return m_list_crack;
  }
  /*}}}*/
  /*observer manuvar{{{*/
  void append(Observer* p_observer) {
    m_list_observer.emplace_back(p_observer);
  }
  void dettach(Observer* p_observer) {
    auto it =
      std::find(m_list_observer.begin(), m_list_observer.end(), p_observer);
    if(it != m_list_observer.end())
      m_list_observer.insert(it, p_observer);
    return;
  }

  void notify() const {
    for(auto p_observer : m_list_observer)
      p_observer->update(this);
    return;
  }
  /*}}}*/
  /*iteration access{{{*/

  typename PointList::iterator getPointBegin() {
    return m_list_point.begin();
  }
  typename PointList::iterator getPointEnd() {
    return m_list_point.end();
  }
  typename PointList::const_iterator getPointBegin() const {
    return m_list_point.begin();
  }
  typename PointList::const_iterator getPointEnd() const {
    return m_list_point.end();
  }

  typename ShapeSpList::iterator getShapeBegin(const size_t& dim) {
    return m_list_shapes.at(dim).begin();
  }
  typename ShapeSpList::iterator getShapeEnd(const size_t& dim) {
    return m_list_shapes.at(dim).end();
  }
  typename ShapeSpList::const_iterator getShapeBegin(const size_t& dim) const {
    return m_list_shapes.at(dim).begin();
  }
  typename ShapeSpList::const_iterator getShapeEnd(const size_t& dim) const {
    return m_list_shapes.at(dim).end();
  }

  typename CrackList::iterator getCrackBegin() {
    return m_list_crack.begin();
  }
  typename CrackList::iterator getCrackEnd() {
    return m_list_crack.end();
  }
  typename CrackList::const_iterator getCrackBegin() const {
    return m_list_crack.begin();
  }
  typename CrackList::const_iterator getCrackEnd() const {
    return m_list_crack.end();
  }
  /*}}}*/
  /*misc {{{*/

  void plot(io::DataPloter&    plt,
            const std::string& seg_cfg = " w l ls 2 lc rgb \"blue\"",
            const std::string& gp_cfg  = " w p pt 0 lc rgb \"black\"",
            const std::string&         = "\0") const {
    m_mesh.plot(plt, seg_cfg, gp_cfg, "\0");
    for(auto& crack : m_list_crack)
      crack.plot(plt, "w lp ls 2 lc rgb \"red\"", "w p pt 0 lc rgb \"cyan\"",
                 "\0");

    if(dimension == 2) {
      plt.appendSetting("set view equal xyz");
      plt.appendSetting("set view 0,0");
    }
    return;
  }
  /*}}}*/

private:
  /**
   * @brief updateBoundary
   *
   * @param boundary
   * @param phase
   */
  template <typename T>
  void updateBoundary(Boundary<T>& boundary, const TimePhase& phase);
  /**
   * @brief build material list from input file
   */
  void buildMaterial() { /*{{{*/
    GEOXFEM_LOG_SEV(info) << "Building materials...";

    m_list_material.clear();

    for(auto& v : m_reader->getMaterials()) {
      ptree                     data = v.second;
      std::string               type = data.get<std::string>("type");
      std::unique_ptr<Material> mat = Material::getFactory().create(type, data);
      size_t                    id  = data.get<size_t>("id");
      mat->setIndex(id);
      GEOXFEM_LOG_SEV(detail) << *mat;
      m_list_material.emplace(id, std::move(mat));
    }
  } /*}}}*/
  /**
   * @brief setup Parameters
   */
  void setupParams() const { /*{{{*/
    if(m_reader->hasValue("calculation.time_step.theta")) {
      _real value =
        m_reader->data().template get<_real>("calculation.time_step.theta");
      time_theta::set_get(value);
    }
    if(m_reader->hasValue("calculation.iteration.tolerance.relative")) {
      _real value = m_reader->data().template get<_real>(
        "calculation.iteration.tolerance.relative");
      relative_tolerance::set_get(value);
    }
    if(m_reader->hasValue("calculation.iteration.tolerance.absolute")) {
      _real value = m_reader->data().template get<_real>(
        "calculation.iteration.tolerance.absolute");
      absolute_tolerance::set_get(value);
    }
    if(m_reader->hasValue("calculation.gravity_scheme.water_density")) {
      _real value = m_reader->data().template get<_real>(
        "calculation.gravity_scheme.water_density");
      water_density::set_get(value);
    }
    if(m_reader->hasValue("calculation.atm.value")) {
      _real value =
        m_reader->data().template get<_real>("calculation.atm.value");
      atm::set_get(value);
    }
    if(m_reader->hasValue("topology.elem_gauss_num")) {
      ptree value = m_reader->data().get_child("topology.elem_gauss_num");
      i_default_gauss_num::set_get(value);
    }
    return;
  } /*}}}*/
  /**
   * @brief buildMesh
   */
  void buildMesh() { /*{{{*/
    GEOXFEM_LOG_SEV(info) << "Building main mesh...";
    geometry::BuildMesh(m_reader->getTopology(), *this, m_mesh);
    return;
  } /*}}}*/
  /**
   * @brief build Time Sequence
   */
  void buildTimeSequence() { /*{{{*/
    GEOXFEM_LOG_SEV(info) << "Building time sequence...";
    m_list_phase.clear();
    _real time = 0;
    for(auto& v : m_reader->getPhases()) {
      const _real during = v.second.template get<_real>("total_time");
      TimePhase   phase(time, during, v.second);
      m_list_phase.emplace_back(phase);
      time += during;
    }
    return;
  } /*}}}*/
  /**
   * @brief buildInitCracks
   */
  void buildInitCracks() { /*{{{*/
    ptree& cracks_data = m_reader->getCrack();
    if(cracks_data.empty())
      return;

    GEOXFEM_LOG_SEV(info) << "Building initial cracks...";
    m_list_crack.clear();
    for(auto& crack_data : cracks_data) {
      Crack* crack = append(Crack{&m_mesh});
      geometry::BuildMesh(crack_data.second, *this, *crack);
    }

  } /*}}}*/

  _real                    m_time;
  size_t                   m_iter;
  PointList                m_list_point;
  std::vector<ShapeSpList> m_list_shapes;
  MaterialList             m_list_material;
  CrackList                m_list_crack;
  DirichletBC              m_dirichlet;
  NeumannBC                m_neumann;
  Mesh                     m_mesh;
  DOFHandler               m_dof_handler;
  std::vector<TimePhase>   m_list_phase;
  size_t                   m_max_pid;
  Reader*                  m_reader;
  ObserverList             m_list_observer;
};  // namespace model

/*solve problem{{{*/
template <typename MethodPolicy,
          typename CouplePolicy,
          typename DimensionType,
          typename EquationPolicy,
          typename SolvePolicy,
          typename LoadPolicy,
          typename degree>
void Model<MethodPolicy,
           CouplePolicy,
           DimensionType,
           EquationPolicy,
           SolvePolicy,
           LoadPolicy,
           degree>::solveProblem() {
  GEOXFEM_LOG_SEV(info) << "Solving problem...";
  TStep last_tStep;
  _real ctime = 0;
  int   nStep = 0;

  for(auto& phase : m_list_phase) {
    GEOXFEM_LOG_SEV(info) << phase << " start...";
    /* 3.1 restore nodal dofs' boundary condition status */
    for(auto vertex : shapesIn(m_mesh, 0))
      for(auto& dof : dofsIn(vertex->toVertex())) {
        dof.second.setBCType(bcType::free);
        dof.second.setTimeFunction(nullptr);
      }
    /* 3.2 update Dofs at boundary condition for current time phase */
    updateBoundary(m_dirichlet, phase);
    updateBoundary(m_neumann, phase);
    /* 3.3 begin time step loop */
    _real etime = phase.getEndTime();
    do {
      m_time      = ctime + phase.getDeltaTime();
      m_time      = (m_time < etime) ? m_time : etime;
      _real dtime = m_time - ctime;
      /* 3.3.1 instanciate current time step */
      TStep curr_tStep(nStep, ctime, dtime);
      solve_time_step<MethodPolicy, CouplePolicy, LoadPolicy,
                      SolvePolicy>::apply(this, curr_tStep);
      /* 3.3.2 output */

      /* 3.3.3 step into next time step */
      last_tStep = curr_tStep;
      ctime      = m_time;
    } while(ctime < etime);
  }
  return;
}
/*}}}*/
/*update boundary information{{{*/
template <typename MethodPolicy,
          typename CouplePolicy,
          typename DimensionType,
          typename EquationPolicy,
          typename SolvePolicy,
          typename LoadPolicy,
          typename degree>
template <typename T>
void Model<MethodPolicy,
           CouplePolicy,
           DimensionType,
           EquationPolicy,
           SolvePolicy,
           LoadPolicy,
           degree>::updateBoundary(Boundary<T>&     boundary,
                                   const TimePhase& phase) {
  GEOXFEM_LOG_SEV(info) << "Updating " << boundary.getClassName() << "...";
  boundary.clear();
  boundary.setPhase(&phase);
  ptree data = boundary.getData();
  for(auto& bc : data) {
    size_t      node_id  = bc.second.get<size_t>("node");
    std::string dof_name = bc.second.get<std::string>("dof");
    Dof*        dof      = getShape(0, node_id)->toVertex().getDof(dof_name);
    if(dof->getBCType() == bcType::dirichlet)
      continue;

    std::string func_name = bc.second.get<std::string>("type") + "_Function";
    auto        t_func    = TimeFunc::getFactory().create(
      func_name, bc.second, phase.getStartTime(), phase.getEndTime());
    dof->setTimeFunction(t_func.get());
    dof->setBCType(T::type);
    boundary.append(dof, t_func);
  }
  return;
}
/*}}}*/

}  // namespace model
}  // namespace geoxfem

#endif /* IMODEL_H */
